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ABSTRACT 

We have included the chemical rate network responsible for the formation of molecu- 
lar Hydrogen in the iV-body hydrodynamic code, Hydra, in order to study the formation 

■ of the first cosmological at redshifts between 10 and 50. We have tested our implemen- 
tation of the chemical and cooling processes by comparing N-body top hat simulations 

■ with theoretical predictions from a semi-analytic model and found them to be in good 
agreement. We find that post-virialization properties are insensitive to the initial abun- 
dance of H 2 . Our main objective was to determine the minimum mass (Msg(z)) of 
perturbations that could become self gravitating (a prerequisite for star formation), 

Q-f and the redshift at which this occurred. We have developed a robust indicator for de- 

tecting the presence of a self-gravitating cloud in our simulations and find that we can 
do so with a baryonic particle mass-resolution of 40M Q . We have performed cosmo- 
logical simulations of primordial objects and find that the object's mass and redshift 
at which they become self gravitating agree well with the Msg(z) results from the top 
hat simulations. Once a critical H 2 fractional abundance of ~ 5 x 10 -4 has formed in 
an object, the cooling time drops below the dynamical time at the centre of the cloud 
and the gas free falls in the dark matter potential wells, becoming self gravitating a 
dynamical time later. 

Subject headings: cosmology: theory — early universe — galaxies: formation — hydro- 
dynamics — molecular processes 



1. Introduction 

In cold dark matter (CDM) cosmologies, the first objects are expected to form between redshifts 
of 10 and 50 (Peebles 1993). Radiative cooling via the rotational and vibrational transitions of 
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molecular hydrogen provides the mechanism by which the baryons may first condense in small dark 
matter potential wells of mass ~ 10 6 M and virial temperature ~ 1000 K; too cool for hydrogen 
line cooling to be effective. Of the several different molecules present in the early universe — for 
example H 2 , HD, and LiH — molecular hydrogen is the most abundant and dominates cooling 
between 100 K and 1000 K. Understanding the cooling and formation of the first objects and their 
influence on subsequent cosmic evolution is an important goal of post-recombination cosmography. 
Of particular interest is the variety of feedback mechanisms which may operate following the epoch 
of first star formation. 

Even a small number of objects at this early time could have a significant effect on the the 
entire universe. A fraction of just ~ 10~ 5 of the baryonic mass condensing into massive stars, 
would produce sufficient energetic photons to re-ionize the universe (Carr et al. 1984; Couchman 
& Rees 1985; Loeb 1997). Observations of the CMB provide a limit on the earliest time that 
reionization could occur. Anisotropics in the CMB on angular scales below the horizon size of 
~ 10° are diminished by free electron production during reionization, and if reionization occurred 
much before z ~ 50, the level of anisotropy in the CMB on ten-degree scales would be lower than 
is observed (Knox 1998; Shaver et al. 1999). 

The mass function of the first stars is very uncertain. Kashlinsky & Rees (1983) argued that the 
main constituents of primordial objects would be low mass stars, with a few very massive objects 
(VMOs) of masses 10 3 — 10 5 M Q forming at later stages. Other authors (Carr et al. 1984; Loeb 
1997) have argued that the lack of cooling by metals increases the Jeans mass and consequently 
star formation is biased towards the production of high mass stars. Some recent numerical studies 
(Omukai & Nishi 1998; Bromm et al. 1999) find that the first generation stars are quite massive. 
Supernovae at the death of high mass stars will contaminate the primordial gas with metals, greatly 
increasing the ability of the gas to cool, as well as injecting large amounts of kinetic and thermal 
energy into the pregalactic medium. 

The H 2 molecule, with a binding energy of only 4.48 eV, is easily dissociated and UV photons 
produced by the first objects could quickly destroy all H 2 , even accounting for self-shielding (Haiman 
et al. 1997). It has been suggested that this might delay further structure formation until larger 
objects form with non-linear mass scales above 10 8 M and virial temperatures of 10 4 K where 
hydrogen line cooling dominates. An alternative view is that the elevated electron abundance 
occurring in the photo-ionized regions around massive stars provides ideal conditions for the re- 
formation of a substantial H 2 fraction as these regions cool and recombine. In the latter case we 
might anticipate that further cooling and condensation would be promoted. 

In order to assess the nature of the first objects and importance of H 2 -mediated cooling, 
many authors have simulated the formation of primordial clouds. The computational expense of 
following serveral chemical rate equations in 3-D hydrodynamic calculations has resulted in the 
use of one-dimensional models of cloud collapse by a number of authors. Lepp & Shull (1984), for 
example, used a spherical model for cloud collapse and followed the chemical reactions responsible 
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for H 2 formation. A similar, high resolution (but still 1-D) calculation was undertaken by Haiman 
et al. (1996). Tegmark et al. (1997) used an even more computationally economical, semi-analytic 
approach which runs so quickly that the collapse mass versus redshift parameter space could be 
extensively explored. 

One-dimensional approaches, however, neglect the complicated filamentary and sheet-like struc- 
tures ubiquitous to three-dimensional simulations of CDM universes. With improvements in nu- 
merical techniques and computer hardware it is now feasible to undertake realistic high resolution 
3-D hydrodynamic simulations which incorporate the relevant chemical reactions. The full three- 
dimensional approach has been taken by Anninos et al. (1997), and is taken here. Despite the recent 
successes of 3-D numerical approaches, significant challenges still remain. Adequate modelling of 
feedback from photio-ionizing radiation and supernova blast waves has not yet been addressed fully 
self-consistently in a numerical simulation although important first steps have been taken (see for 
example Gnedin & Ostriker 1997 and Haiman et al. 1999). Second, the shallow slope of the CDM 
spectrum on scales < 10 8 M leads to a wide range of scales collapsing nearly co-asvally and very 
long-range correlations in the peaks of the mass fluctuation field, posing a severe challenge to the 
limited dynamic range of 3-D simulations. In this paper we present numerical techniques for mod- 
elling the collapse and condensation of gas in dark matter haloes at high redshift with a view to 
reliably identifying the sites of first star formation. In a subsequent paper we will investigate the 
effect of feedback of these first stars on the pregalactic medium. 

A pre-requisite for star formation in dark-matter-dominated cosmologies is that gas becomes 
self-gravitating (i.e. the baryon density becomes locally greater than the dark matter density) 
although this is perhaps not a sufficient condition. The earliest time that stars could then form is 
the time at which there exist objects that have self-gravitating regions; likely in the core. Identifying 
these regions in numerical simulations is the key objective of this paper. Our goal is to model 
primordial objects and determine their mass and the time at which they become self-gravitating 
within cosmological simulations. This aim distinguishes our approach from that of Bromm et al. 
(1999) in which the focus is to examine the internal properties of the first self-gravitating clouds. 

Structure formation is hierarchical in the CDM model: small objects collapse first and coalesce 
into larger objects. We may indentify three main stages in the merger hierarchy of the two- 
component fluid as follows (cf. White & Rees 1978). Initially, the potential gradients resulting 
from the first small dark matter fluctuations will be very much less than pressure gradients which 
can be generated in the gas, leading only to small adiabatic density perturbations in the baryonic 
component. As the dark matter haloes become larger through merging and accretion, the potential 
gradients will increase to a point at which the baryonic material is compressed to the same fractional 
overdensity as the dark matter. A dark matter halo is conventionally said to have reached the Jeans 
mass when this occurs. As the mass of the underlying dark matter haloes increases beyond this 
point the gas becomes subject to the dominant gravitational influence of the virializing dark matter 
haloes. This will result in the baryonic matter shock-heating to the effective virial temperature of 
the dark matter halo and attaining a broadly similar density profile. The subsequent behaviour of 
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gas in a dark matter halo depends upon the efficiency with which it can cool. In hierarchical models 
the ratio of the cooling time to the dynamical time is of crucial importance (White & Rees 1978; 
Blumenthal et al. 1984). If t coo \ <C td yn , the gas cools efficiently and free falls in the dark matter 
potential well. Conversely, if t coo \ 3> td yn , the gas is unable to condense appreciably before being 
incorporated into a larger halo as the hierarchy builds, a process which, in a smooth hierarchy, 
occurs on a characteristic timescale set by td yn ■ In this paper we demonstrate explicitly these three 
stages as haloes grow in the CDM hierarchy. 

Note that the Jeans mass as used in the above sense for a two component fluid differs from 
its original definition as the perturbative scale separating acoustic oscillation from instability in a 
self-gravitating gas. A perturbative treatment of the two component fluid is complicated in detail 
(see, for example, Meiksin et al. 1999), but we may make a useful simple estimate of the mass 
scale of a virialized (non-linear) dark matter halo below which baryonic matter is not significantly 
perturbed, by equating the velocity dispersion of the dark matter in the virialized halo with the 
specific thermal energy of the gas: 

1 , 2 \ 3GM 3kT 

2>> = ^iT = W (1) 

At the Jeans mass the gas has been adiabatically compressed to the same overdensity as the dark 
matter, implying a baryon temperature T = Tb(5 v + l) 2 / 3 , where Tj, is the background baryon 
temperature and 5 V is the overdensity relative to the background at the time of collapse. The Jeans 
mass is then approximately 



firriH J 2GHq 
where Tq(1 + z) 2 is the temperature of the unperturbed gas at redshift z. 

The outline of the paper is as follows. We discuss the network of chemical reactions involved 
in the production of molecular hydrogen in § 2, along with the details of its implementation in 
our iV-body gravitational and gasdynamic code, Hydra (Couchman et al. 1995). As an initial test 
of our implementation we compare, in § 3, the results of spherically symmetric top hat collapses 
with an analytic model. In § 4 we determine the minimum mass at a given redshift that an object 
must have to cool efficiently and collapse to high densities where star formation may begin. This 
was accomplished by performing many different simulations of top hat perturbations with different 
masses and overdensities. A more realistic, cosmological simulation is described in section § 5. 
Here we simulate a (25kpc) 3 volume of space and follow the behaviour of the baryonic component 
in detail through the three stages of cooling and condensation in hierarchical models as described 
above. 
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2. Molecular hydrogen cooling and chemistry 

The H 2 molecule has long been known to be an important coolant for gravitational collapse 
of hydrogen clouds that are cooler than T ~ 1000 K. Gould & Salpeter (1963) discussed forma- 
tion processes and cooling rates for molecular hydrogen in the interstellar medium; Gould (1964) 
discussed protostellar collapse through molecular hydrogen cooling. H 2 is also a key ingredient 
of the globular cluster formation scenario of Peebles & Dicke (1968). The hydrogen molecule is 
also important for structure formation in the early universe; the first objects to form have virial 
temperatures of ~ 1000 K, where the H 2 molecule is the dominant coolant. Through H 2 cooling, 
star clusters (Couchman Sz Rees 1985) or VMOs (Bond et al. 1984) could form and subsequently 
reionize the intergalactic medium. Molecular hydrogen is also essential to the formation of the first 
generation stars (Stahler 1986). 

In addition to hydrogen, the primordial gas contained Li, D, 3 He, and 4 He with abundances of 
10~ 10 , 10~ 5 , 10~ 5 , and 0.24 respectively, relative to hydrogen. Line cooling is dominated by neutral 
hydrogen because the excitation rates of Li, D, and 3 He are much smaller (Abel et al. 1997). 

Molecular hydrogen is the most common molecule at the epoch at which the first primordial 
stars form (Lepp & Shull 1984). In the ISM molecular hydrogen formation occurs on the surface 
of grains, but, since stars produce the grains, this cannot be an important mechanism for the first 
objects. There are, however, other chemical reactions that are able to form H 2 , which are discussed 
below. At a redshift of 100, the fractional abundance of molecular hydrogen was 10 -6 (Galli & 
Palla 1998). Other molecules that were important coolants are HD and LiH which had abundances 
10~ 3 and 10~ 4 times lower than that of H 2 . The HD molecule becomes the dominant coolant at 
20 K < T < 100 K, and LiH dominates below 20 K. These two molecules will be important in 
the coldest and densest gas in an object, and to get to that state the gas initially will radiate 
thermal energy through H 2 cooling. Molecular hydrogen cooling is efficient enough to produce self 
gravitating objects, which is a necessary condition for star formation. Since we are interested in 
following the evolution of primordial objects only to the self gravitating state where pb « pd m an d 
not to higher baryonic densities, we need include only H 2 cooling and can ignore HD and LiH 
cooling. 

2.1. Cooling processes 

When an H 2 molecule collides (usually with neutral hydrogen, the most abundant species) 
and becomes rotationally or vibrationally excited, it may radiatively de-excite which results in 
gas cooling, or it may collisionally de-excite which results in no net energy loss from the gas. 
H 2 vibrations are excited at 6000 K; below 3000 K only the rotational levels of the molecules can 
be excited (Puy et al. 1993). Radiative de-excitation dominates at the low densities relevant for 
primordial gas. In this case, H 2 molecules are mainly in the ground or first (J = 1) rotational 
state. Collisional excitations are quickly followed by radiative decay. At high densities (n = 10 3 
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to 10 4 cm 3 for 10 2 K <T< 10 4 K), collisions dominate and Boltzmann populations are reached. 
These densities are achieved in the late stages of the collapse of primordial objects. 

The cooling mechanisms at various temperatures are shown in figure 1. At temperatures 
between 10 2 K and 10 4 K the H 2 molecule is the most effective coolant. Hollenbach & McKee 
(1979) express the H 2 cooling function Ah 2 in the form 

A H2 [n H ,T] = ^ LTE — , (3) 
1 + n cr /n H 

where Alte is the LTE cooling function (see Hollenbach h McKee 1979 for the full expression), 
critical density defined by: 

n cr _ h-LTE 

n H A H2 [n H -> 0] ' 

and Ah 2 [riff — ► 0] is the low-density limit of the cooling function. The H 2 cooling function has 
been computed by Martin et al. (1996) for T > 600 K and Forrey et al. (1997) for T < 600 K. Galli 
& Palla (1998) provided a fit to these two computations in the low-density limit, which we adopt 
here: 

log A H2 [n H -> 0] = -103.0 + 97.59 log T - 48.05(logT) 2 + 10.80(logT) 3 - 0.9032(log T) 4 . (5) 



There is uncertainty in the H 2 cooling function because of the difficulty in calculating the 
interaction potential at low temperatures. Different choices for rotational and vibrational H-H 2 
rate coefficients will produce differences in Ah 2 ; figure 2 illustrates H 2 cooling rates computed by 
several authors. 

At temperatures between 10 4 K and 5 x 10 4 K radiational de-excitation of H dominates the 
cooling function. We use the expression given by Dalgarno & McCray (1972): 

A Hline ~ 7.5 x 10- 19 e - 118348K / T n 2 / c -(l - / e -). (6) 

Cooling processes such as He line cooling and Bremsstrahlung (fig. 1) that operate at higher tem- 
peratures may be ignored, since the primordial objects have virial temperatures of ~ 1000 K. 



2.2. Molecular hydrogen formation and destruction 

Despite a relatively small number of elements in the primordial gas, there is a complicated 
network of reactions involving these elements; Galli & Palla (1998) list 87 reactions in their com- 
prehensive study of the chemistry of the early universe. Fortunately many of these do not affect 
significantly the collapse of primordial objects and the number of important reactions is much 
smaller. The HD and LiH molecules are not important at densities in our regime of interest. He- 
lium chemistry and cooling becomes important above 10 5 K, which is well above the 1000 K virial 
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Fig. 1. — Cooling rates for Bremsstrahlung (dotted line), H (dashed line) and He (dash-dotted line) 
line cooling, and H 2 (circles) cooling. The e~, H + , He + , and He ++ abundances were computed 
assuming collisional equilibrium, and the H 2 fractional abundance was 3 x 1CP 4 , which is typical 
for early objects. 
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Fig. 2. — H 2 cooling rates per H 2 molecule in the low density (n < 0.1 cm -3 ) limit, from: Galli 
k Palla 1998 (solid line), Martin et al. 1996 (stars), Tegmark et al. 1997 (circles)., Hollenbach & 
McKee 1979 (dash-dotted line), and Lepp & Shull 1984 (dashed line). 
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Table 1. Reactions involved in the formation and destruction of molecular hydrogen. 



Reaction Rate Reference 

(1) H+ + e" — ► H + 7 1.88 x 10- 10 r-°- 64 2,1(11) 

(2) H + e" — ► H + + 2e" see reference 1(11) 

(3) H + e" — ► H-+7 1.83 x 10- 18 T - 88 2,1(12) 

(4) H-+H — ► H 2 + e" 1.3 xlO" 9 3,1(10) 

(5) H-+e" — ► H + 2e" 1.6 x 10~ 6 (1 + 34/T 5 / 6 ) exp(-65/T 1 / 3 ) 3,1(10) 

(6) H 2 + e" — ► 2H + e" 5.6 x lO^T 1 / 2 exp(-102124/T) 1(13) 

(7) H 2 + H+ — ► H+ + H 2.6 x 10- 8 exp(-2.12/T 4 )/T 1 / 2 8,1 

(8) H-+H+ — ► 2H 4xlO" 6 / Tl/2 7,1(15) 

(9) H + H — ► H 2 4 x 10~ 27 @ T = 100 K 3 

(10) H+ + H — ► H 2 + H+ 6.4 xlO" 10 1(5) 

(11) H+ + H- — ► H 2 + H < 2.3 x 10- 7 3,1(15) 

(12) H+ + e~ — ► H 2 + H 5 x 10~ 9 6 

(13) H 2 + H 2 — ► H 2 + 2H 3 x 10~ 4 exp(-5.2/T 4 )/T 3 / 2 3 

(14) H 2 + H~ — ► H 2 + H + e" very small 3 

(15) H 2 + H+ — > H+ + H 2.1 x 10~ 9 3 

(16) H 2 + e" — > H + H" 2.7 x 10- 8 exp(-4.3/T 4 )/T 3 / 2 3 

(17) H 2 + H — > 3H see reference 1(14) 

(18) H+ + e" — ► 2H 4.2 x 10- 8 /TV 2 6,1(16) 

(19) H+ + H — > H+ + 7 1.85 x 10- 23 T 18 1(4) 

(20) H~ + H+ — > H+ + e" 2.291 x lO" - 4 for T < 2 x 10 4 K 1(4) 

(21) H~+H — > 2H + e" see reference 1(10) 

(22) H + 7 — ► H+ + e" see reference 1(17) 

(23) H-+7 — > H + e" 0.114T 2 - 13 exp(-8650/T 7 ) 9,1(4) 

(24) H+ + 7 — ► H + H+ 6.36 x 10 5 exp(-71600/T 7 ) 9,1(18) 

(25) H^ + 7 — > 2H+ + e _ see reference 1(4) 

(26) H 2 + 7 — ► + e_ see reference 1 

(27) H 2 + 7 — > H+ + H + e~ see reference 19 

(28) H 2 + 7 — ► H(ls) +H(2s,2p) see reference 1 

(29) H 2 + 7 — > 2H(ls) see reference 1 



References. - 1: Abel et al. 1997, 2: Hutchins 1976, 3: Hirasawa 1969, 4: Shapiro & Kang 
1987, 5: Karpas et al. 1979, 6: Solomon & Werner 1971, 7: Hill & Silk 1975, 8: Rapp & Francis 
1962, 9: Tegmark et al. 1997, 10: Janev et al. 1987, 11: Ferland et al. 1992, 12: Wishart 1979, 
13: Donahue & Shull 1991, 14: Dove & Mandy 1986, 15: Dalgarno & Lepp 1987, 16: Schneider 
et al. 1994, 17: Osterbrock 1974, 18: Standi 1994, 19: Couchman 1986 

Note. — References in parenthesis are the source used by Abel et al. (1997). Where more than 
two references are given, the reaction rate listed was obtained from the first. 
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temperatures that primordial objects have. Ignoring HD, LiH, and He molecular species reduces 
the number of relevant reactions considerably, and will only limit investigations of the latest stages 
of baryonic condensation at low temperatures which are not the focus of this paper. We summarize 
and justify in this section the reactions which it is necessary to include in our numerical model to 
accurately follow the chemistry and cooling of the first objects. 

Molecular hydrogen production proceeds via only two major mechanisms. In the channel 
the proton acts as a catalyst: 

H+ + H — ► H++7 
H+ + H — ► H 2 + H+. 

In the H~ channel the electron acts as a catalyst: 

H + e" — ► H" + 7 
H~ + H — ► H 2 + e". 

For redshifts higher than 200, CMB photons are energetic enough that the reaction H~ + 7 — ► 
H + e~ neutralizes all H ions and therefore the H + channel dominates for z > 200. At lower 
redshifts, the H~ channel proceeds much more quickly than the H + channel and dominates H 2 
production. 

The destruction of the intermediaries H - and H;j~ proceeds much more quickly than their 
production. This leads to equilibrium fractional abundances of H~ and (relative to total 
hydrogen number density) of: 

, = fafnfc- + k i6fu 2 f c - / 7 -j 

H hfn + k 8 fu+ + hf c - + hif u + + hofm + foifu 

, k 7 fu 2 fu+ + fcl9/H+ + fe20/H~/H+ (8] 

fclO + fcll/H- + We- ' 

where the rate for reaction i is fcj. Table 1 lists the reactions we have considered here; where more 
than one source is listed for the reaction rate, we have used the rate given by Abel et al. (1998). 

Reaction rates depend on temperature and the abundances of the reactants, all of which change 
with time, and therefore these must be known to identify the dominant reactions. This leads to a 
circular argument: the density and temperature must be known to identify dominant reactions, but 
the dominant reactions must be known to properly model the density and temperature evolution. 
One could include every chemical reaction, but this is computationally expensive. Instead, we used 
reasonable estimates for the chemical abundances, density, and temperature to tentatively identify 
the important reactions, and included these in our TV-body code, along with the dominant cooling 
mechanisms; the implementation is discussed in §2.5 below. Then we performed a spherically 
symmetric top hat collapse of 10 6 M test object (see §3 for details) following the temperature and 
abundances during the evolution to plot the reaction rates versus redshift in figures 3-5. From these 
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figures the unimportant reactions may be identified, and their exclusion justified. We believe that 
our identifications are robust. 

Figure 3 shows the H 2 production and destruction mechanisms along with the temperature 
of the collapsing object as a function of redshift. The peak temperature occurs at 1 + z = 20 
when the test object virializes. Although this redshift was arbitrarily chosen, the behaviour of the 
reactions does not change significantly for perturbations that virialize between 50 < z < 10, the 
redshifts at which the first objects are expected to form in models with CDM-like spectra. It is 
immediately obvious from figure 3 that the only H 2 production mechanism of importance at these 
redshifts is the H~ channel. The abundance is less than 10~ 16 (Galli & Palla 1998) so reaction 
12 is unimportant. Since the abundance does not influence the H 2 production significantly, 
reactions 10, 11, and 18 may be ignored. 

Figure 4 shows the reactions involving the H~ ion. Reactions 3 and 4 (the H~ channel) are 
in equilibrium and no other reactions proceed at comparable rates. At T > 10 4 K the H~ ion is 
quickly destroyed by collisions with electrons (reaction 5), thus terminating H 2 production through 
this route. 

At the low temperatures (T < 5000 K) encountered in the the first bound objects, there is 
no efficient means of destroying H 2 . Although unimportant for objects with virial temperatures of 
T < 1000 K, H 2 dissociation mechanisms must be included for the reaction network to be accurate 
at higher temperatures. The five H 2 destruction mechanisms in figure 3 all have similar rates; 
however at T = 5000 K protons begin to dissociate H 2 via reaction 7. Once the free electron 
abundance becomes significant at T > 2 x 10 4 K, H 2 dissociation by electrons dominates, with the 
production of two neutral hydrogen atoms and an electron (reaction 6) strongly favoured over the 
production of a neutral hydrogen and a negative hydrogen ion (reaction 16). The dissociation of H 2 
by H (reaction 17) is unimportant since almost all hydrogen is ionized at T > 2 x 10 4 K, at which 
point H 2 is dissociated entirely by collisions with free electrons (reaction 6). At lower temperatures 
dissociation of H 2 by H is unimportant since charge exchange with H + (reaction 7) occurs at a 
greater rate. H 2 -H 2 and Hg-H^ collisions are rare because of their low abundance so reactions 13 
and 15 are negligible. 

Since H 2 production depends on the free electron abundance, reactions involving the ionization 
of hydrogen must also be examined. The neutralization of H + by an electron (reaction 1) is in 
equilibrium with ionization of H by an electron (reaction 2) in figure 5, and these two reactions 
control the degree of hydrogen ionization. 

The dominant reactions thus 1 through 7. All reactions involving helium have been disregarded, 
which is permissible so long as the temperature is restricted to T < 3 x 10 4 K where helium begins 
ionization and becomes an important coolant. All molecules other than H 2 (such as LiH and HD 
and those involving He) have not been included since their abundances are negligible. 

The equations that describe the evolution of the various species can be determined using the 
dominant reactions discussed above. The total particle abundance is n = n[H] + n[H + ] + 2n[H 2 ] 
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Fig. 3. — Reaction rates for the formation and destruction of molecular hydrogen. Temperatures 
and fractional abundances of reactants (/a and /#) used to compute the rates were taken from a 
top hat simulation of an early object collapsing at 1+z = 20. H 2 creation via the H~ channel clearly 
dominates. The five H 2 destruction reactions all progress at similar rates at these temperatures. 
The numbers in parentheses indicate the reaction in table 1. 
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Fig. 4. — Reaction rates for the formation and destruction of H~. The two reactions in the 
H channel for molecular hydrogen production (top line) are in equilibrium. At T > 10 4 K the 
neutralization of H by e~ becomes important. These three reactions exclusively control the H~ 
abundance. The numbers in parentheses indicate the reaction in table 1. 
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Fig. 5. — Reaction rates involving ionized hydrogen. Radiative recombination controls the H + 
abundance. The numbers in parentheses indicate the reaction in table 1. 
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and the fractional abundance of species A is f A = n[A]/n. The H, H + , and H 2 fractional abundances 
evolve as: 

= k ifn+f c - ~ hfnfe- (9) 



n dt 
1 dfm 
n dt 
1 dfu 2 
n dt 



= k 2 fuf e --k 1 f H+ f e - (10) 
= kifn-fn - k 7 fu 2 fn+ - hfH 2 f e - (11) 



Since the H formation and destruction reactions are in equilibrium it is never necessary to integrate 
the chemical equations to determine its abundance. Note that the H~ equilibrium abundance 
(equation 7) may be simplified by including only the rate coefficients for the dominant reactions 
(i.e. k 3 , fe 4 , fc 5 ). 



2.3. Photoionization and Photodissociation 

Although ionization and dissociation by photons are important in the early universe, by z = 100 
there are essentially no CMB photons with energies above the thresholds for reactions 22 (13.6 eV) 
and 23 (0.755 eV). Since becomes important in molecular hydrogen production only at redshifts 
z > 200, reactions 24 and 25 may be neglected. The dissociation reactions 26, 27, 28, and 29 are 
inhibited because the CMB photons are much less energetic than the H 2 binding energy of 4.48 eV. 
Since there are no sources of energetic photons, all photoionization and photodissociation reactions 
may be ignored from z = 100 until after the first objects have been formed. 



2.4. Shocks 

There are several issues concerning H 2 chemistry that must be addressed when the gas shocks 
during collapse. The H 2 molecule has a small binding energy of 4.48 eV, and if a shock has a 
velocity of v > ^Jl x 4.48 eV/(2m p ) ~ 20km/s then H 2 may be collisionally dissociated (Kwan 
1977). The first objects to form will have masses between approximately 10 5 M Q and 10 6 M 
and virial temperatures of T ~ 1000 K. This temperature corresponds to an average thermal 
velocity v = \J (3kT/2)/(m/2) which is 5km/s for H and 3.5km/s for H 2 , well below the collisional 
dissociation velocity. We have not observed any cases in which the shock velocity approaches the 
dissociation velocity in our simulations of primordial objects. 



2.5. Incorporation of the chemical reactions into Hydra 



Our iV-body gravitational and gasdynamic code, Hydra, is an adaptive mesh implementation 
of the smoothed particle hydrodynamics (SPH) algorithm and the particle-particle particle-mesh 
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(P 3 M) algorithm. A complete description of the code is given by Couchman et al. (1995). The 
chemical state of the gas is carried by assigning each SPH particle a fractional abundance of H 2 
and H + as described below. 

If the macroscopic changes caused by gravitational and adiabatic heating or cooling operate 
much more slowly than the chemical processes, then the two become decoupled and may be inte- 
grated separately. This is indeed the case for our simulations. During a spherical top hat collapse, 
the maximum values of T/T and p/p are on the order of 10 7 and 10 8 years, while /h+//h+ is 
approximately 10 4 years. 

The time rate-of-change of density is computed for each particle via p = (pi — pi-\)/dt where 
pi and are the density values at the current and previous time step. The value of p is assumed 
to be constant over the time step; while this is not strictly true, the approximation is justified 
since the time steps are typically a factor of 100 smaller than the characteristic time scale p/p. 
The time rate-of-change of temperature from adiabatic processes is computed analogously: T a( i^ = 
(T a d,i — Tad^i-i) I 'dt. These two slowly varying dynamical quantities, p and T a( i, are needed to 
compute the rapidly- varying chemical evolution of the system. The cooling functions (fig. 1), 
adiabatic processes and the chemical abundance equations (9, 10 and 11) are then simultaneously 
integrated for each particle at every time step. We assume that the H + and e~ abundances are 
equal, which is justified provided temperatures are constrained to T < 10 5 K. The H~ abundance 
may be computed using equation (7). The abundance is never needed since it is not a reactant in 
any of the dominant reactions. The neutral hydrogen abundance is computed via the conservation 
equation n = n[H] + n[H + ] + 2n[H 2 ]. Of the six species present in the dominant reaction network 
(e~, H, H~, H + , H 2 , H^), only two are independent, and each particle therefore carries values for 
the H 2 and H + abundances. 

We use the stiff ODE integrator STIFBS from Numerical Recipes (Press et al. 1992). This 
routine was also chosen for a very similar problem by Haiman et al. (1996) who experimented with 
other Numerical Recipes routines and found STIFBS to be the most efficient. They also tested the 
LSODAR routine of Hindmarsh (1983) and found STIFBS ~ 20% faster while producing identical 
results. We compared our integration of the chemical network of ODEs to a stiff ODE integrator in 
the Matlab software package and found them to agree to 1 part in 10 4 , which is more than adequate 
for our purposes. 
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3. Tests of the code 

3.1. Comparison with analytic models 

In this section we present two tests: i) top hat evolution from quiet initial conditions which 
proceeds to a singular collapse; and ii) a top hat collapse in which perturbations within the top hat 
overdensity lead to virialization at roughly half the maximum (turnaround) radius. 

The top hat model has two free parameters, the initial mass Mth an d the initial overdensity 
which determines the redshift at which the analytic top hat collapses to a singularity. Using an 
ODE integrator in the Matlab software package, we numerically integrated the top hat density 
evolution along with relevant chemical reactions and cooling functions to analytically determine 
the chemical and temperature evolution of spherical perturbations. 

For a simple initial test of the chemical network in Hydra we simulate a spherically symmetric 
perturbation of uniform density and compare that to the analytic top hat evolution. Particles arc 
placed into the simulation box on a regular lattice, and a spherical region enclosing Mth is cut 
out and compressed by an amount that causes the singularity to be reached at the desired redshift. 
Although placing particles on a regular lattice is unphysical (because it leads to a singularity), 
the absence of noise in the particle distribution allows an exact comparison between the code and 
the analytic solution. The particles within the top hat are assigned velocities according to the 
time derivative of the cycloid equations. The gas particles are given an increase in temperature 
corresponding to the adiabatic compression. 

Figure 6 shows the dark matter, baryon, temperature, H~, and H 2 evolution computed using 
the analytic top hat model and our iV-body code for a top hat perturbation of total mass 1O 6 M . 
The simulation was terminated prior to z = 20 since there is no point in following the evolution 
through the singularity. This is however more than sufficient to test the accuracy of the chemistry 
in the iV-body code. We performed many other simulations using various masses and virializa- 
tion redshifts, and there was always excellent agreement with the semi- analytic model up to the 
singularity. 

In physical reality the singularity is avoided since random perturbations in the density field will 
cause non-radial motions to develop, and instead of collapsing to a point the object will virialize. 
The analytic evolution for the top hat density and temperature therefore must be modified to 
approximate the virial behaviour. Tegmark et al. (1997) modeled the formation of primordial 
objects with a semi-analytic approach. They assumed that the density followed the top hat solution 
until the temperature reached the virial value, after which the density remained constant at 187T 2 
times the background density at the nominal collapse redshift. If the overdensity reached the virial 
value before T > T v i r , they held the density constant and assumed that the gas would be shock 
heated to the virial temperature. This density function, the cooling rate, and chemical reactions 
were integrated to follow the temperature and chemical abundances in the top hat. 
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Fig. 6. — Comparison of temperatures, baryonic and dark matter densities, H + , and H 2 fractions 
given by an analytic top hat model (solid lines) and our iV-body code, Hydra (symbols) for a 
perturbation of total mass 10 6 M that virializes at 1 + z = 20. 
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Many cosmological simulations use a regular cubic lattice for the initial particle distribution, 
onto which is imposed a power spectrum of fluctuations according to the cosmology. For a simple 
top hat model, the only fluctuation imposed is a simple compression of a spherical region of space. 
A cubic grid is inappropriate in this case since collisionless particles would collapse to a singularity 
and oscillate about the collapse centre. We introduce noise into the particle distribution in order 
to generate substructure which will lead to virialization. 

To make a fair comparison between the semi-analytic top hat model and the simulations, the 
initial conditions must be set up to mimic the semi-analytic behaviour as closely as possible. This 
however involves two antagonistic criteria. Very smooth initial conditions will not produce as much 
substructure and will follow the semi- analytic top hat model closely up to the virialization redshift. 
However, the smooth initial conditions do not allow the dark matter to virialize well: much of the 
dark matter rebounds outward to a large distance while also producing a very dense, small virialized 
region which has an overdensity well above the semi-analytic value of 5 ~ 200. Adding more noise 
to the initial conditions promotes the formation of substructure which results in a larger portion 
of the dark matter contained in the virialized region and a smaller central overdensity, but the 
increased substructure causes the temperature and density to rise slightly above the semi-analytic 
value prior to virialization. 

The amount of substructure that forms before the dark matter virializes is governed by the 
noise in the initial particle positions and the time elapsed between the start of the simulation and 
virialization. Without altering the initial particle positions, substructure in the dark matter may 
be quantitatively increased by simply starting the simulation at an earlier redshift. One must now 
decide upon an appropriate amount of substructure to use. We have chosen to use random initial 
conditions and start the simulations at a redshift that is 20 expansion times prior to virialization (eg. 
Zi = 400 for z v = 20; these values are appropriate for a top hat containing roughly 10 4 particles). 
This arrangement produces a virialized region with approximately 80% of the dark matter within 
the analytic virial radius, and a mean overdensity within the half-mass radius of approximately 
200. This creates post-virialization conditions that are similar to the semi-analytic model, at the 
expense of pre-virialization agreement. Baryons will fall into the shallow potential wells provided 
by the substructure present in the dark matter, slightly elevating the density, temperature, and 
H 2 abundance. The amount of disagreement with the semi analytic model prior to virialization 
is small and has negligible impact upon the post-virialization behaviour. Requiring the presence 
of substructure is also attractive since this reflects more closely collapses in a hierarchical CDM 
cosmology — they are never monolithic, as they are in the analytic top hat model. 

We performed several top hat runs and compared with the semi-analytic predictions. We find 
that after virialization, for objects that remain pressure supported, the simulations agree well with 
the semi-analytic model because the baryon overdensity remains near the virial value. For objects 
that are able to cool and become self gravitating, the baryon overdensities become much higher than 
the virial value, and the semi-analytic model does not agree well with the simulations. However, if 
the post-virialization baryon density in the semi-analytic model was raised to the value observed 
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in the simulations, the chemical abundances and temperature agreed very well. 



3.2. Sensitivity to initial chemical abundances 

The H 2 abundance in the early Universe has been computed by many authors and varies 
considerably: « 1CT 4 (Tegmark et al. 1997; Haiman et al. 1996), 5 x 10~ 7 (Lepp & Shull 1983), 
7.5 x 10~ 6 (Puy et al. 1993). The simulations performed here use the H 2 abundance given by 
Anninos & Norman (1996): 

3 /2 

/H 2 =2xlO" 2 °M^(l + z ) 5 - 1 (12) 

where zq is the largest redshift at which H 2 can be efficiently formed through the channel. 
Peebles (1993) found the residual ionization after recombination to be: 

/ H + = 1.2 x l(T 5 ft 1/2 /(M^)- (13) 



To determine the sensitivity that the initial H 2 abundance has on the collapse, two iV-body 
simulations of top hat perturbations were run with initial H 2 abundances of 2 x 1CT 6 and ~ 10 -4 . 
Figure 7 shows the evolution of the baryonic and dark matter densities, the temperature, and 
the H 2 and H + abundances for the central region of a top hat with M = 10 6 M that virializes at 
1 + z v = 20. The H 2 abundance for z < z v is not sufficient in either run to alter the gas temperature 
and hence the gas evolves adiabatically: it cools with the expansion up to the turnaround redshift 
zt = 30 and subsequently heats as the perturbation begins to contract. After turnaround the 
increase in baryonic density accelerates the molecular hydrogen production and the difference in 
the H 2 abundances between the two runs quickly disappears. Both runs reach the same final value 
°f /h 2 = 10~ 3 shortly after virialization. Another run was done using an initial H 2 abundance of 
10~ 8 which shortly stabilized at /h 2 = 5 x 10~ 6 . Prior to virialization, there are no reactions that 
efficiently destroy H 2 , and the factor limiting the H 2 abundance is the number of free electrons 
available to catalyze the H~ channel. The temperature, density, and chemical evolution after 
virialization are thus insensitive to the initial molecular hydrogen abundance. 

The molecular hydrogen abundance reaches a plateau before virialization of 6 x 1CT 6 which 
disagrees with the simulations of Haiman et al. (1996) and Tegmark et al. (1997) who found the 
plateau values to be 1.6 x 10~ 4 and 10 -4 respectively. These authors have used a rate for the 
photodissociation of H;j~ (important at z > 200) that is too low, which consequently caused a 
higher H 2 abundance, as noted by Abel et al. (1998). Our H 2 abundance agrees with the newer 
calculations of Galli & Palla (1998) and Abel et al. (1998). 

Regardless of the redshift at which the simulations are started, the chemical abundances are 
initialized at z = 100 and from then on the chemical reactions and cooling processes are followed. 
Computing the abundance of H 2 at earlier redshifts necessitates inclusion of more chemical reaction 
which would slow the code, only to produce an H 2 abundance at z = 100 which is already known 
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(see e.g. Galli & Palla 1998). At this redshift, the top hat simulations are still quite smooth and the 
baryonic density and temperature are not significantly different from the background values. Some 
error in the initial abundances is undoubtedly committed by neglecting to compute the chemistry 
at redshifts earlier than 100, but as shown above, variations in the initial H 2 abundance do not 
alter the post-virialization behaviour. 

3.3. Softening issues 

In order to avoid excessively large forces between two particles in close proximity, the gravi- 
tational attraction is "softened" so that the gravitational force decreases at small distances. The 
distance at which the softened gravity diverges from the Newtonian value is set by the softening 
length, s. This parameter must be carefully tuned to the problem, otherwise the results will be 
tainted by unphysical numerical effects. If s is too small then two body collisions become appar- 
ent; if s is too large then the softening may prevent a group of particles from reaching a density 
high enough to correctly model a self-gravitating region — a necessary condition for star formation. 
SPH forces are computed by using a smoothing kernel to average over the approximately 32 nearest 
neighbours; the SPH smoothing length (h) is chosen adaptively such that it encloses that number of 
neighbours. In very dense regions, however, it is necessary to constrain the SPH smoothing length 
to a minimum value (one half of the gravitational softening length) to balance the gravitational 
and hydrodynamic resolutions, otherwise unphysical clustering of the baryon particles can occur 
(Thacker et al. 1998). The minimum smoothing length determines the maximum density that the 
gas particles may reach. Efficient cooling leads to dense clustering of the baryons relative to the 
dark matter particles, and a mismatch in the ideal s and h m i n . We have performed several tests to 
quantify the effects that different softening and smoothing lengths produce. 

The overdensity of an object in which gas pressure is initially unimportant is, at virialization, 
1 + 5 = 187r 2 in an Q = 1 universe (Gunn & Gott 1972). The mean separation between particles 
in the perturbation is then d = (lSvr 2 ) -1 / 3 ^ where df, g is the mean interparticle separation of the 
uniform background matter. For a simulation box volume of V = L 3 containing N = n 3 particles, 
dbg = L/n and d ~ 0.2L/n. The softening length should enclose at least roughly 10 particles 
to minimize spurious two body effects, which gives s = 2d ~ OAL/n. This may be taken as an 
estimate of the minimum softening length, since the central region of a virialized object has an 
overdensity well in excess of <5 ~ 200. As an extreme, if S = 10 5 , then s ~ 0.04L/n. 

The smoothing length has a much greater effect on the baryonic density than on the dark 
matter density for cases in which the baryons can cool efficiently. This is illustrated in figure 8 
which shows the baryonic and dark matter density profiles for top hat simulations using various 
smoothing lengths (table 2, runs 1 to 4). In these runs, enough H 2 is formed that the gas can cool 
efficiently and free fall into the dark matter potential well. If s is too large, then the baryonic density 
Pb remains lower than the dark matter density pdm- The smoothing length has the undesirable 
ability to determine if the gas can become self gravitating (i.e. pb > Pdm)- This problem does not 
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Fig. 7. — A comparison of two runs having initial molecular hydrogen abundances of 2 x 10 -6 
(thick lines) and ~ 10 (thin lines). The final physical and chemical states of the object after 
virialization (z = 20) are insensitive to the initial H 2 abundance. 
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exist for hot, pressure supported objects — in these, even an extremely small softening will not 
cause the object to become self-gravitating (figure 9). 

Since we are concerned with finding the minimum mass of an object that can cool efficiently 
and become self-gravitating, the situation above where the softening has such a strong influence 
on the baryonic density is unsatisfactory. A collapse criterion sensitive to numerical (not physical) 
parameters is not robust. In a pressure-supported object, the number of particles within the 
softening radius is approximately 32, while in a self gravitating object the number of particles will 
be magnitudes larger. Figures 10 and 11 show the mass and temperature profiles for an object 
that has cooled and an object that is pressure supported. The number of particles (and also mass 
and temperature) within the softening radius of an object that has cooled stays fairly constant 
regardless of the size of the softening radius. As an alternative to the comparison between pb and 
Pdm comparison, we consider an object to be have collapsed to a self gravitating state if there 
are more than one thousand particles within the softening radius at a temperature of just over 
100 K. In practice it is quite easy to distinguish between objects that have become self gravitating 
and those that are pressure supported since the former display rapid and obvious changes in their 
density and temperature profiles. 



3.4. Particle resolution 

Increasing the number of particles in an iV-body simulation provides better resolution, at 
the expense of computational speed. Since we intended to perform many simulations of top hat 
perturbations of different masses and virialization redshifts in order to determine the mass required 
for an object to become self gravitating at a certain redshift, it was desirable to find the minimum 
number of particles that would provide acceptable resolution. This minimum number will also 
provide a pointer to the particle mass resolution necessary for robust simulations of the first objects 
in a cosmological setting. 

We performed several simulations using various particle resolutions, and found that the central 
density in N = 2 x 16 3 runs was lower than that in N = 2 x 32 3 runs. Some variation with resolution 
has to be expected since the central regions in the lower resolution runs are marginally resolved. 
In objects where the central gas cools quickly, within the softening radius there are about 50 
dark matter and 500 gas particles in the 2 x 16 3 simulations. The 2 x 32 3 runs have well defined 
cores containing approximately 600 dark matter and 6000 gas particles within the softening radius. 
Throughout the paper we will use the term "core" to describe the regions in self gravitating objects 
that have constant (with respect to distance to centre) high density and low temperature. 

The N = 2 x 32 3 runs in which the gas is able to cool efficiently have N gas (r) (cumulative 
number of gas particles) profiles that are flat, which indicates that almost all of the gas particles 
initially contained in the top hat perturbation lie within the softening radius. In this case the 
resolution is limited by the softening length and not the number of particles, but for the N = 2 x 16 3 
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Fig. 8. — Density profiles of an object that can cool efficiently (runs 1 to 4) for various softening 
lengths. Shown are the dark matter density (thin lines — computed dividing the mass in a radial 
shell by its volume) , and the baryon density (thick lines — computed by averaging the SPH density 
of all particles within a radial shell). The vertical lines show the softening lengths used for each 
simulation. The runs with the two smallest softening values have self-gravitating cores (p& > pdm)- 
Increasing the softening further causes the baryonic density to drop below that of the dark matter 
density. 
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Fig. 9. — Density profiles of a pressure supported object (runs 5 to 8) for various softening lengths. 
In the pressure supported case, softening does not determine if the central region becomes self- 
gravitating (pb > pdm)- As in figure 8, thick lines are baryon densities and the thin lines are dark 
matter densities. The vertical lines show the softening lengths used for each simulation. Densities 
beyond the softening length are convergent with densities computed using a smaller softening length. 
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Fig. 10. — Radial profiles of gas temperature (solid lines) and gas mass (broken lines) for various 
smoothing lengths (runs 1 to 4). The object contains about 9000 particles, has a total mass of 
10 6 M.q and is able to cool efficiently. As the smoothing length (vertical lines) is increased the peak 
temperature decreases, but the temperature of the baryons in the core of the object (r < lpc) is 
not significantly influenced by the smoothing length. There are approximately 3000 gas particles 
within the smoothing length regardless of its value, and they have cooled to a temperature where 
H 2 is no longer an effective coolant. 
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Fig. 11. — Radial profiles of gas temperature (solid lines) and gas mass (broken lines) for various 
smoothing lengths (runs 5 to 8) for an object of total mass of 2 x 10 5 M which is not able to cool 
efficiently. As the smoothing length is decreased (vertical lines), the central temperature decreases 
and approaches the central temperature observed in an object that is able to cool efficiently (fig. 
10), however the central density is very much less for pressure supported objects. 
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run the situation is reversed. This demonstrates that the minimum limit of acceptable particle 
resolution has been reached. We performed one top hat simulation with N = 2 x 64 3 and found the 
results agreed with the N = 2 x 32 3 run and are confident therefore that 2 x 32 3 particles provide 
adequate resolution. 

The dark matter alone drives the gas dynamics until the baryonic material becomes self grav- 
itating. The differences in the distribution in the baryonic matter between the N = 2 x 16 3 and 
N = 2 x 32 3 simulations is primarily caused by the differences in the dark matter distribution. 
The mass profiles in the N = 2 x 32 3 runs were close to M{r) oc r 3 as expected for a perturbation 
of uniform density, but in the N = 2 x 16 3 runs there simply are not enough particles to accu- 
rately model the smooth density of the top hat model, and the mass profiles deviate from being 
proportional to volume. 
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Table 2. Simulation parameters for softening and resolution tests 



run 


mass 


box size 


softening 




N 




(10 6 M ) 


(kpc) 


(pc) 






1 


1.03 


20.8 


3.5 


20 


2 x 32 3 


2 


1.03 


20.8 


7.0 


20 


2 x 32 3 


3 


1.03 


20.8 


14.0 


20 


2 x 32 3 


4 


1.03 


20.8 


28.0 


20 


2 x 32 3 


5 


0.17 


11.4 


1.9 


20 


2 x 32 3 


6 


0.17 


11.4 


3.7 


20 


2 x 32 3 


7 


0.17 


11.4 


7.5 


20 


2 x 32 3 


8 


0.17 


11.4 


14.7 


20 


2 x 32 3 
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4. Determination of minimum mass threshold for baryonic condensation 

In CDM-like cosmologies, smaller dark matter perturbations collapse at higher redshifts, but 
for the baryonic matter to become self gravitating, enough H 2 must exist for the gas to be able to 
cool efficiently. The rate of H 2 formation is related to the baryonic density which is initially driven 
by the dark matter potential well. A minimum mass exists above which the H 2 production rate 
is sufficient to cool the gas efficiently. This mass threshold is a function of redshift, and in this 
section we explore this parameter space by simulating top hat perturbations of various masses. Our 
main interest is to perform more realistic cosmological simulations of primordial object formation, 
and the analytic and top hat simulations provide a reasonable estimate of the mass scale at which 
collapse should occur at a given redshift. All simulations presented here are in an Einstein-de Sitter 
universe with f^m = 0.94, Of, = 0.06, and a Hubble constant of Hq = GSkms™ 1 Mpc -1 . 

An object must contain a certain mass of dark matter to be able to force the baryons into 
a self gravitating state where star formation may begin. This mass is redshift dependent: as the 
universe expands it becomes less dense, and deeper dark matter gravitational potential wells are 
required to drive up the baryonic density to a point where H 2 production becomes efficient enough 
that the object may cool on a dynamical time. In this section, using both the semi-analytic and 
iV-body methods, we have determined the minimum mass required for the baryons in an object to 
become self gravitating (Msg) as a function of redshift. 

Figure 12 shows the minimum mass threshold versus redshift relation originally computed by 
Tegmark et al. (1997). Using the same method, our result differs considerably since we have used 
a different molecular hydrogen cooling function (see fig. 2), and have included different chemical 
reactions. Tegmark et al. (1997) used the cooling function of Hollenbach & McKee (1979), and 
assumed that the H 2 molecules are all in the para state and that H 2 vibrational cooling was 
negligible. At temperatures higher than 1000 K, H 2 vibrations may be excited and vibrational 
cooling becomes important; this is why the cooling function that Tegmark et al. (1997) used drops 
below the other curves when T > 1000 K. More efficient cooling allows smaller, less dense objects 
to collapse, and consequently our computation of Msg(z) is a factor of approximately 3 smaller 
at z = 50. At z = 10, Tegmark et al. (1997) found that the objects able to cool efficiently had 
temperatures of approximately 10 4 K, and at this temperature their cooling function is two orders 
of magnitude lower than the cooling function we used. Also, at temperatures greater than 10 4 K 
H 2 dissociation reactions become important. Since they did not incorporate any H 2 destruction 
mechanisms their H 2 abundances are too high, though this is not enough to offset the less efficient 
cooling function. We find that objects able to cool efficiently at z = 10 are 25 times less massive. 

The results of several N = 2 x 32 3 top hat simulations are also shown for comparison with 
the analytic model (fig. 12). Inside the top hat there are ~ 2 x 10 4 particles, and in cases where 
gas in the central region cools quickly, the condensed cores are well resolved and contain ~4x 10 3 
particles. For a top hat of total initial mass Mth destined to virialize at z = z v , we have plotted 
a star if the collapse succeeded (using the collapse criterion developed in section 3.3) within one 
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dynamical time after virialization and a circle if the object remained pressure supported. Some 
disagreement between the simulations and the semi-analytic result is to be expected since the 
latter assumes that the baryonic overdensity remains constant at 5 = 200 after virialization, and 
is uniform within the top hat. In spite of the assumptions made in the semi-analytic result and 
the different criteria used to determine if the collapse was successful, the two methods agree to 
within a factor of 2. The slope of the Msg(z) line determined with the semi-analytic method is 
slightly steeper than for the simulations but we do not consider this difference to be of special 
significance. Ultimately the precise determination of Msg{z) is dependent upon the details of the 
chemical reaction rates and cooling rates, and the amount of substructure present in the object. 
Further, although the threshold for collapse is very sharp because of the density-cooling feedback 
(a small increase in mass can cause an appreciable increase in density), the precise position of the 
line separating collapse from pressure support is highly sensitive to small variations in reaction 
rates and small numerical changes and hence is simply indicative of the general behaviour in the 
M-z plane. Our results show that we expect the collapse mass to be robust within a factor of 2 at 
a given redshift. 
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Fig. 12. — Minimum mass threshold required for self gravitation (Msg) versus redshift. The dash- 
dotted line shows the result of Tegmark et al. (1997), and the dashed line shows our semi-analytic 
result. The dotted line is the Jeans mass. The symbols show results of iV-body simulations using 
2 x 32 3 particles. Objects less massive than Msg at a given redshift will remain pressure supported 
(circles); more massive objects cool efficiently and the baryonic material in the central region 
becomes self-gravitating (stars). The solid line is an approximate fit to the simulation results. 
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5. Cosmological simulation of first objects 

Our objective here is to follow the growth of objects through merging and accretion in the 
CDM hierarchy and to be able to identify objects that are able to form enough H 2 that they cool 
efficiently and become self gravitating. Of primary interest is the determination of the time at 
which the object exceeds the Jeans mass threshold, and the time at which the cooling time drops 
below the dynamical time. To follow objects through these phases we simulated a (25/i _1 kpc) 3 
volume in an 17 = 3.4, f^/fj = 0.06, h = 0.65 universe. An overdense universe was used simply 
to model the environment of a 3 — 4u perturbation which will lead to the formation of the first 
cosmological objects. The simulation contained equal numbers of gas and dark matter particles, in 
total N = 2 x 64 3 . This gives a mass of 5.1 M© for the baryon particles and 81 M for the dark 
matter particles. 

The mass evolution and merger history of the most massive object in the simulation is shown 
in figure 13. We define M200 to be the mass enclosed by a thin shell that has an overdensity of 
200. The complete family tree is very large and to reduce the complexity of the plot, at each 
redshift only the direct ancestors of the most massive object is shown. Lines connect the smaller 
objects that have merged into the largest object. The separation between the self-gravitating and 
pressure-supported domains is shown for both the semi-analytic method and the N = 2 x 32 3 top 
hat simulations. The agreement between the cosmological simulation and the top hat simulation 
computation of Msg is quite good; the point at which the object becomes self gravitating is close 
to the top hat result. 

In Figure 14 we plot the entropy evolution of the approximately 2000 gas particles which reside 
in the core of the final cooled object. This figure shows clearly the three stages in hierarchical 
evolution. Before the dark matter halo reaches the Jeans mass, the gas will respond adiabatically 
to the growing dark matter potential well, and thus the entropy (defined as oc T/p 2 / 3 ) will remain 
constant at the background value. Once the Jeans mass is reached infalling gas will be shock heated 
(a non-adiabatic process) and hence will have its entropy raised. Once cooling becomes effective, 
the entropy in the core drops rapidly. 

The baryonic matter in objects less massive than the Jeans mass is pressure supported against 
gravitational contraction. In these objects the baryonic matter should be relatively unperturbed by 
the dark matter, and will not have steep baryonic radial density profiles. As the baryonic material 
in an object approaches the Jeans mass, the baryonic radial density profile will begin to steepen 
and become similar in shape to the dark matter radial density profile. The most massive object in 
the simulation exhibits this behaviour, as shown in figure 15. By z = 28 the baryon density profile 
has started to steepen but its slope is still less than that of the dark matter. At z = 24 the object 
has exceeded the Jeans mass and the slopes of the baryon and dark matter density profiles become 
very close. 

At z = 28, the temperature and H 2 fraction (relative to the total hydrogen abundance) in the 
most massive object produce a cooling time that is significantly greater than the dynamical time 
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Fig. 13. — The mass evolution and merger history are shown for the most massive object in the 
cosmological simulation. The thin lines connect objects to their predecessors. The dotted line 
shows Msg computed with the semi-analytic method, and the solid line shows the results from the 
N = 2 x 32 3 top hat simulations. Initially, the baryonic matter in the objects is pressure supported 
(circles), and as mass is accreted the gas becomes self gravitating (stars). Physical properties of 
the indicated objects are shown in figures 15, 16, and 17. Within r2oo there are 5100, 2500, and 
280 dark matter particles and 4600, 2100, and 160 gas particles at redshifts of 20, 24, and 28 
respectively. 
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Fig. 14. — The evolution of the entropy of the ~ 2000 core particles in the most massive object 
in a cosmological simulation. The solid line shows the mean entropy of the core particles, dotted 
lines show the particles with the minimum and maximum entropy, and the dashed line shows one 
standard deviation above the mean. Initially the gas is subject to adiabatic processes only, and the 
entropy does not change with time. As the gas begins to fall into the dark matter potential well 
and exceeds the Jeans mass, it is shock heated and the entropy begins to rise. After enough H 2 
has formed, the cooling time falls below the dynamical time, rapidly reducing the entropy, and the 
gas in the core becomes self gravitating within a dynamical time. 
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Fig. 15. — The evolution of the radial profiles of the mass and overdensity for baryonic (solid lines) 
and dark matter (broken lines) are shown for the most massive object in a cosmological simulation 
at redshifts of 28, 24, and 20. These profiles correspond to the object at the epochs indicated in 
figure 13. At z = 24 enough H 2 has formed that the baryons cool efficiently and free fall in the 
dark matter potential well. By z = 20 — 10 7 years or approximately one dynamical time later - 
the object becomes self gravitating. 
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Fig. 16. — The cooling time (dash-dotted lines), dynamical time (dashed lines), and molecular 
hydrogen fraction (solid lines) are shown at redshifts of 28 (thin lines) and 24 (thick lines) for the 
most massive object in a cosmological simulation. From z = 28 to z = 24 the increase in density 
and temperature accelerate the H 2 production. The increase in temperature and H 2 abundance 
conspire to reduce the cooling time below the dynamical time in part of the central region of the 
object by z = 24. 
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Fig. 17. — The evolution of the temperature profile of the baryons in the most massive object in 
a cosmological simulation at redshifts of 28, 24, and 20. These profiles correspond to the objects 
indicated in figure 13. The central temperature decreases as time evolves. Shock heating of gas 
falling into the deepening dark matter potential causes the peak temperature to increase with time, 
and, as the accretion shock moves outward, the peak temperature occurs at greater distances from 
the centre of the object. 
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(fig. 16). Up to this time, the gas has been relatively mildly perturbed by the growing dark matter 
halo. As the dark matter mass in the object increases, the deepening gravitational potential well 
drives up the gas density and consequently accelerates H 2 production. By z = 24, the H 2 abundance 
and temperature have increased enough that near the centre of the halo the cooling time drops 
below the dynamical time which is 10 7 years at that time. The cooling time is greater than the 
dynamical time very close to the centre of the object, however here the gas pressure and temperature 
are limited by the softening length, which is approximately 2 pc at this time. The physical state 
of the object is such that collapse is inevitable: the dark matter provides a deep enough potential 
well, and there is enough molecular hydrogen that the gas cools efficiently and free falls in the 
dark matter potential well. This increases the gas density and the molecular hydrogen production, 
thereby accelerating the H 2 cooling, and a runaway feedback mechanism begins. One dynamical 
time later, the core temperature is just above 100 K (fig. 17), and the baryonic density exceeds the 
dark matter density and the object becomes self gravitating (fig. 15). At this point the conditions 
in the cold, dense environment are appropriate for star formation. 

To test for possible numerical resolution effects, we re-sampled the initial conditions of the 
N = 2 x 64 3 simulation to produce an N = 2 x 32 3 distribution of particles with similar random 
density fluctuations. The radial profiles of various physical quantities for the most massive object 
were compared for both runs and found to agree well. This convergence indicates that 2 x 32 3 
particles provides sufficient numerical resolution, at least for the most massive object. This test 
suggests that a mass resolution of 40 M Q per baryonic particle and 650 M Q per dark matter particle 
is sufficient to adequately identify the first self gravitating clouds in a cosmological simulation. 

Although we have discussed only one object here, we have simulated several realizations of 
the same cosmological model to ensure that the object presented was not anomalous in any way. 
Since the two different resolutions discussed above were found to be convergent, we used the lower 
resolution and varied only the random number seed to give different particle positions. The few 
most massive objects in these simulations became self gravitating at redshifts ranging from 30 to 15, 
and all exhibited behaviour similar to that described above: a few dynamical times after exceeding 
the Jeans mass, t coo \ fell below td yn and the objects became self gravitating within approximately 
one dynamical time. Although the redshifts at which the objects became self gravitating covers a 
large range, the objects had different masses and all agreed with the M$g{ z ) prediction of the top 
hat simulations to within a factor of two. 

Since there exists some uncertainty in the H 2 cooling rate and chemical reaction rates, we 
investigated the effect on the collapse of primordial objects of varying key rates. Recent deter- 
minations of the H 2 cooling function have differences of a factor ~ 2, arising from uncertainty in 
the rotational and vibrational H-H 2 rate coefficients (Galli & Palla 1998). Various authors also use 
different rates for the formation of molecular hydrogen. We performed two simulations that differed 
only in the H 2 cooling rate and the rate of H formation (which governs the total H 2 abundance) . 
Inflating Ah 2 and the H~ formation rate each by a factor of two caused the most massive object to 
become self gravitating only slightly earlier (z = 19 versus z = 17) and the overall impact on the 
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evolution of the object was quite small. 

By comparing the H 2 cooling time scale to the dynamical time, Tegmark et al. (1997) argued 
that objects will collapse if virial temperatures are high enough to produce a critical H 2 fraction of 
~ 5 x 10~ 4 . Abel et al. (1998) recovered the same abundance in self gravitating regions of objects 
that they simulated. However, they argued that the H 2 fraction is dependent upon Ah 2 , and that 
the use of different cooling functions will produce different values for the H 2 abundance threshold 
needed for collapse. They stated therefore that their agreement in the critical H 2 fraction was 
coincidental since they used different molecular hydrogen cooling functions. We have also found 
approximately the same H 2 fraction present in the cores of self gravitating objects (fig. 16). The 
cooling function that we have used (Galli & Palla 1998) is quite similar to that used by Tegmark 
et al. (1997) for T < 400 K, so the agreement between the H 2 abundances is not surprising. The 
agreement between our work and that of Abel et al. (1998) does seem unexpected, since the cooling 
function differs markedly from the one we employed. We offer a possible explanation. A more 
efficient cooling function will help to maintain the number of free electrons to catalyze the H~ 
channel, but will also retard the production of H~ since the rate falls with temperature. These two 
effects conspire to minimize differences in the H 2 caused by variations in the H 2 cooling function. 
To investigate this quantitatively, we performed used several different cooling functions in the semi- 
analytic method and recovered each author's results. Thus, although different cooling functions 
will produce different temperatures, the resulting change in the H 2 abundance is small. 
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6. Summary 

We have studied the general framework of cooling and condensation of primordial objects in 
two component hierarchical cosmologies. Molecular hydrogen is the dominant coolant in primordial 
objects with virial temperatures between 100 K and 1000 K. Incorporating molecular hydrogen 
cooling into a numerical model requires explicit integration of a chemical network because of the 
non-equilibrium abundances of several reactants. The relevant H 2 formation and cooling processes 
have been added to the cosmological gravitational and SPH code, Hydra. We have focused on the 
question of identifying haloes within which baryons can cool efficiently and in which the first stars 
will form, but do not attempt to model the internal structure of these first clouds in detail. 

Tegmark et al. (1997) used a simple semi-analytic approach to model the density of a top hat 
perturbation, and integrated this along with the chemical equations and H 2 cooling function to 
determine the evolution of the perturbation. This procedure has been followed here with a recently 
computed H 2 cooling function as a test of the TV-body code. The semi-analytic model approximates 
virialization by assuming that the perturbation reaches an overdensity of 200, remaining at constant 
physical density thereafter. We find excellent agreement between the numerical method and the 
semi-analytic model for cloud collapses in which the gas remains pressure supported. In cases in 
which efficient cooling occurs, the baryon density may substantially exceed the naive analytic value 
and the numerical values diverge. If, however, the baryon overdensity at virialization in the semi- 
analytic model is adjusted to match the N-body results, the species abundances come into good 
agreement once again. We find that the final, post-virialization, species abundances, densities and 
temperatures are insensitive to the high redshift abundance of H 2 . 

In accordance with the results of Tegmark et al. (1997) we find a sharp threshold (denoted by 
M$g) i n the dark matter halo mass-redshift plane separating haloes in which gas can cool from 
those in which it cannot. However, our use of a current cooling function yields Msg(z) values up to 
an order of magnitude smaller than found in the earlier work. The threshold found by the iV-body 
method and the semi-analytic method is found to agree to within a factor of 2 over the redshift 
range considered. 

Our primary purpose is to identify haloes in which the gas may become self gravitating. For 
this to occur, the baryons have to condense by a linear factor ~ (Sl/f^) 1 / 3 , leading to a disparity in 
the ideal softening lengths for the dark matter and gas. In order to maintain similar gravitational 
and hydrodynamic resolution, and to avoid spurious numerical effects, we do not allow the softening 
or hydrodynamic resolution to follow arbitrarily high baryon densities but instead rely on counting 
the number of SPH particles within an SPH smoothing length. We find this technique to be a 
robust indicator of whether a halo contains self gravitating gas or not. By this means we are able 
to adequately detect the formation of the first self gravitating gas in CDM-like cosmologies with a 
gas particle mass of ~ 40 M . 

We have performed cosmological simulations of primordial object formation in a CDM universe. 
The evolution and merger histories of several objects were traced through the simulation and we 
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were able to show that the descriptive picture of baryon condensation in dark haloes proposed 
by White & Rees (1978) is an accurate representation of the relevant processes leading to the 
condensation of gas in dark haloes. We have shown that H 2 cooling is very efficient: once the H 2 
abundance reaches a level such that the cooling time in the centre of a halo is less than the local 
dynamical time, the gas free falls in the potential well becoming self gravitating one dynamical 
time later. The point at which the cooling time first drops below the dynamic time agrees very well 
with the minimum mass threshold given by the top hat simulations. The critical H 2 abundance 
required for collapse, 5 x 10 -4 , is the same as found by Tegmark et al. (1997). We find that this 
abundance is insensitive to the details of the cooling function used. 
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